// Author: Ivan Chau
package uk.ac.cam.kcic2.math;

import java.awt.geom.QuadCurve2D;
import java.util.ArrayList;

public class Polynomial {

	private final static int QUADRATIC_ROOTS = 2;

	private final double[] coeffs;
	private ArrayList<Double> realSolutions = new ArrayList<Double>();

	public Polynomial(double[] coefficients){
		coeffs = new double[coefficients.length];
		for (int i = 0; i < coefficients.length; i++){
			coeffs[i] = coefficients[i];
		}
		switch(coeffs.length){
		case 3:	realSolutions = quadraticSolver();	break;
		//case 4: realSolutions = cubicSolver();	break;
		//case 5: realSolutions = quarticSolver();	break;
		default: throw new IndexOutOfBoundsException();
		}
	}

	public double[] getRoots(){
		double[] roots = new double[realSolutions.size()];
		if (realSolutions.size() > 0){			
			for (int i = 0; i < realSolutions.size(); i++){
				roots[i] = realSolutions.get(i);
			}			
		}
		return roots;
	}

	private ArrayList<Double> quadraticSolver(){
		// coeffs = {c,b,a} for ax^2 + b^x + c = 0		
		final double a = coeffs[2];
		final double b = coeffs[1];
		final double c = coeffs[0];
		ArrayList<Double> realRoots = new ArrayList<Double>();

		// separate the case when b and/or c = 0 because solveQuadratic
		// does not return the correct no. of roots which are zeros
		// check for c = 0
		if (Math.abs(c) < Params.EPILSON){
			realRoots.add(0d);
			if (Math.abs(b) < Params.EPILSON){
				realRoots.add(0d);
				return realRoots;
			}
			else{
				realRoots.add(-b/a);
				return realRoots;
			}
		}
		else{
			double[] roots = new double[QUADRATIC_ROOTS];
			int numOfRoots = QuadCurve2D.solveQuadratic(coeffs,roots);
			if (numOfRoots >= 1){
				realRoots.add(roots[0]);
				if (numOfRoots == 2){
					realRoots.add(roots[1]);
				}
			}
			return realRoots;
		}
	}
/*
	private ArrayList<Double> cubicSolver(){
		// coeffs = {d,c,b,a} for ax^3 + bx^2 + c^x + d = 0
		ArrayList<Double> realRoots = new ArrayList<Double>();
		final double a = coeffs[3];
		final double b = coeffs[2];
		final double c = coeffs[1];
		final double d = coeffs[0];

		// check for the trivial case of a quadratic
		if (Math.abs(a - 0d) < Params.EPILSON){
			return quadraticSolver();
		}

		// otherwise, depress the cubic:
		// sub x = t - b/(3a) and divide equation by a to give
		// t^3 + pt + q = 0
		final double p,q;
		final double[] t = new double[CUBIC_ROOTS];
		if (!(Math.abs(a - 1d) < Params.EPILSON && Math.abs(b - 0d) < Params.EPILSON)){
			p = (3d*a*c - b*b) / (3d*a*a);
			q = (2d*b*b*b - 9d*a*b*c + 27d*a*a*d) / (27d*a*a*a);
		}
		else{
			p = c;
			q = d;
		}

		// triple real root - t^3 + q = 0
		if (Math.abs(p - 0d) < Params.EPILSON){
			double xRoot = Math.cbrt(-q) - b/(3d*a);	// use x = t - b/(3a)
			for (int i = 0; i < CUBIC_ROOTS; i++){
				realRoots.add(xRoot);
			}
			return realRoots;
		}
		// three real roots
		if (4d*p*p*p + 27d*q*q <= 0d){
			// solve using trig functions if all roots are real
			// sub t = ucos(theta) in t^3 + pt + q = 0
			// choose u such that 
			// we have 4cos3(theta) - 3cos(theta) - cos(3theta) = 0
			// choose u = 2*sqrt(-p/3) and divide equation by u^3/4
			// gives 4cos3(theta) - 3cos(theta) - 3q/(2p) sqrt(-3/p) = 0
			// so we have cos(3theta) = 3q/(2p) sqrt(-3/p)
			// hence, t = u*cos((1/3)*arccos(3q/(2p) sqrt(-3/p) - k*2pi/3)), k=0,1,2
			// let trigSolveCubicT(p,q) be t[0]
			// and using -pi <= arccos(u) <= pi, -1 <= u <= 1
			t[0] = trigSolveCubicT(p,q);
			t[2] = -trigSolveCubicT(p,-q);
			t[1] = -t[0] - t[2];
			// use x = t - b/(3a)
			for (int i = 0; i < CUBIC_ROOTS; i++){
				realRoots.add(t[i] - b/(3d*a));
			}
			return realRoots;
		}
		// one real root and a pair of complex conjugate roots
		// two possible situations:
		// solve using hyperbolic functions
		if (4d*p*p*p + 27d*q*q > 0d && p < 0d){
			// solve using cosh
			t[0] = -2d*Math.abs(q)/q*Math.sqrt(-p/3d)*
					Math.cosh(1d/3d*
							arccosh(-3d*Math.abs(q)/(2d*p)*
									Math.sqrt(-3d/p)));
			// use x = t - b/(3a)
			realRoots.add(t[0] - b/(3d*a));
			return realRoots;			
		}
		if (p > 0d){
			// solve using sinh
			t[0] = -2d*Math.sqrt(p/3d)*Math.sinh(1d/3d*
					arcsinh(3d*q/(2d*p)*Math.sqrt(3d/p)));
			// use x = t - b/(3a)
			realRoots.add(t[0] - b/(3d*a));
			return realRoots;
		}
		else{
			// this line should never be reached,
			// all cases should have been handled explicitly
			throw new NumberFormatException();
		}
	}

	private double trigSolveCubicT(double p, double q){
		double t = 2*Math.sqrt(-p/3d)*Math.cos(1d/3d*
				Math.acos(3d*q/(2d*p)*
						Math.sqrt(-3d/p)));
		return t;
	}

	private double arccosh(double value){
		double x = value;
		if (Math.abs(x - 1) < Params.EPILSON){
			x = 1d;
		}
		if (x >= 1d){
			return Math.log(x + Math.sqrt(x*x - 1d));
		}
		else{
			throw new NumberFormatException();
		}
	}

	private double arcsinh(double x){
		return Math.log(x + Math.hypot(x,1d));
	}

	private ArrayList<Double> quarticSolver(){
		// coeffs = {e,d,c,b,a} for ax^4 + bx^3 + cx^2 + dx + e = 0
		ArrayList<Double> realRoots = new ArrayList<Double>();
		final double a = coeffs[4];
		final double b = coeffs[3];
		final double c = coeffs[2];
		final double d = coeffs[1];
		final double e = coeffs[0];

		// check for cubic
		if (Math.abs(a - 0d) < Params.EPILSON){
			return cubicSolver();
		}

		// otherwise, the job is done if we can factor the quartic into
		// the product of two quadratics, so
		// x^4 + bx^3 + cx^2 + dx + e = 0 becomes
		// (x^2 + px + q)(x^2 + rx + s) = 0 which expands to give
		// x^4 + (p+r)x^3 + (q+s+pr)x^2 + (ps+qr)x + qs = 0
		// thus, b = p+r, c = q+s+pr, d = ps+qr, e = qs
		// and for depressed quartic where b = 0,
		final double c2;
		final double d2;
		final double e2;
		if (!(a == 1d && b == 0d)){
			// sub x = x - b/(4a) and divide the equation by a
			double temp_c2 = -3d*b*b/(8d*a*a) + c/a;
			c2 = Math.abs(temp_c2 - 0d) < Params.EPILSON? 0d : temp_c2;
			double temp_d2 = Math.pow(b/(2d*a),3d) - b*c/(2d*a*a) + d/a;
			d2 = Math.abs(temp_d2 - 0d) < Params.EPILSON? 0d : temp_d2;
			double temp_e2 = -3d/256d*Math.pow(b/a, 4d) + c*b*b/(16d*a*a*a) - 
					b*d/(4d*a*a) + e/a;
			e2 = Math.abs(temp_e2 - 0d) < Params.EPILSON? 0d : temp_e2;
		}
		else{
			c2 = c;
			d2 = d;
			e2 = e;
		}
		// resulting in x^4 + cx^2 + dx + e = 0
		// and so, r = -p, c+p^2 = s+q, d/p = s-q, e = qs
		// eliminate s, q by
		// (c+p^2)^2 - (d/p)^2 = (s+q)^2 - (s-q)^2 gives
		// (c+p^2)^2 - (d/p)^2 = 4e
		// if we set P = p^2, then this becomes
		// P^3 + 2cP^2 + (c^2 - 4e)P - d^2 = 0
		// and solve this cubic to find p
		final double cubicA = 1d;
		final double cubicB = 2d*c2;
		double temp_cubicC = c2*c2 - 4d*e2;
		final double cubicC = Math.abs(temp_cubicC - 0d) < Params.EPILSON? 0d : temp_cubicC;
		final double cubicD = -d2*d2;
		Polynomial bigP = new Polynomial(new double[]{cubicD, cubicC, cubicB, cubicA});
		double[] bigPRoots = bigP.getRoots();
		int numOfBigPRoots = bigPRoots.length;
		ArrayList<Double> depressedRoots;
		switch (numOfBigPRoots){
		case 0: return realRoots;
		case 1:
			if (bigPRoots[0] >= 0d){
				depressedRoots = findRootsFromP(Math.sqrt(bigPRoots[0]), c2, d2);
				if (depressedRoots.size() > 0){
					for (int i = 0; i < depressedRoots.size(); i++){
						// use x = x - b/(4a)
						realRoots.add(depressedRoots.get(i) - b/(4d*a));
					}
				}

			}
			return realRoots;
		case 3:
			for (int i = 0; i < numOfBigPRoots; i++){
				if (bigPRoots[i] >= 0d){
					depressedRoots = findRootsFromP(Math.sqrt(bigPRoots[i]), c2, d2);
					if (depressedRoots.size() > 0){
						for (int j = 0; j < depressedRoots.size(); j++){
							// use x = x - b/(4a)
							realRoots.add(depressedRoots.get(j) - b/(4d*a));
						}
					}
					return realRoots;
				}				
			}
			return realRoots;

		default: throw new NumberFormatException();
		}
	}

	private static ArrayList<Double> findRootsFromP(double p, double c, double d){
		double r,s,q;
		if (Math.abs(p - 0d) < Params.EPILSON){
			r = -p;
			s = c/2d;
			q = c/2d;
		}
		else{
			r = -p;
			s = (c + p*p + d/p) / 2d;
			q = (c + p*p - d/p) / 2d;
		}
		// the quartic x^4 + bx^3/a + cx^2/a + d^x/a + e/a = 0 can be factored as
		// (x + px + q)(x + rx + s) = 0
		Polynomial quadratic1 = new Polynomial(new double[]{q,p,1d});
		Polynomial quadratic2 = new Polynomial(new double[]{s,r,1d});
		ArrayList<Double> roots = new ArrayList<Double>();
		for (double root : quadratic1.getRoots()){
			roots.add(root);
		}
		for (double root : quadratic2.getRoots()){
			roots.add(root);
		}
		return roots;
	}*/
}
